rgdal : interface entre R et les librairies GDAL (Geospatial Data Abstraction Library) et PROJ4.
sp : classes et methodes pour les données spatiales dans R.
rgeos : accès à la librairie d’opérations spatiales GEOS (Geometry Engine - Open Source) : area, perimeter, distances, dissolve, buffer, overlap, union, contains…
Les fonctionnalités de sp, rgeos et rgdal dans un package unique.
Manipulation plus aisée, objets plus simples
Auteur principal et maintainer : Edzer Pebesma (auteur de sp)
Compatible avec les syntaxes pipe et les opérateurs du tidyverse.
Format des objets spatiaux sf
Import de données
library(sf)
mtq <- st_read("data/mtq/martinique.shp")
Reading layer `martinique' from data source `/home/tim/Documents/prz/flo/data/mtq/martinique.shp' using driver `ESRI Shapefile'
Simple feature collection with 34 features and 23 fields
geometry type: POLYGON
dimension: XY
bbox: xmin: 690574.4 ymin: 1592426 xmax: 736126.5 ymax: 1645660
epsg (SRID): 32620
proj4string: +proj=utm +zone=20 +datum=WGS84 +units=m +no_defs
Affichage de données
plot(st_geometry(mtq))
Extraire les centroïdes
mtq_c <- st_centroid(mtq)
plot(st_geometry(mtq))
plot(st_geometry(mtq_c), add=TRUE, cex=1.2, col="red", pch=20)
Construire une matrice de distances
mat <- st_distance(x=mtq_c,y=mtq_c)
mat[1:5,1:5]
Units: m
[,1] [,2] [,3] [,4] [,5]
[1,] 0.000 35297.56 3091.501 12131.617 17136.310
[2,] 35297.557 0.00 38332.602 25518.913 18605.249
[3,] 3091.501 38332.60 0.000 15094.702 20226.198
[4,] 12131.617 25518.91 15094.702 0.000 7177.011
[5,] 17136.310 18605.25 20226.198 7177.011 0.000
Agréger des polygones
mtq_u <- st_union(mtq)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2, border = "red")
Construire une zone tampon
mtq_b <- st_buffer(x = mtq_u, dist = 5000)
plot(st_geometry(mtq), col="lightblue")
plot(st_geometry(mtq_u), add=T, lwd=2)
plot(st_geometry(mtq_b), add=T, lwd=2, border = "red")
Réaliser une intersection
m <- rbind(c(700015,1624212), c(700015,1641586), c(719127,1641586),
c(719127,1624212), c(700015,1624212))
p <- st_sf(st_sfc(st_polygon(list(m))), crs = st_crs(mtq))
plot(st_geometry(mtq))
plot(p, border="red", lwd=2, add=T)
mtq_z <- st_intersection(x = mtq, y = p)
plot(st_geometry(mtq))
plot(st_geometry(mtq_z), col="red", border="green", add=T)
Construire des polygones de Voronoi
google: “st_voronoi R sf” (https://github.com/r-spatial/sf/issues/474 & https://stackoverflow.com/questions/45719790/create-voronoi-polygon-with-simple-feature-in-r)
mtq_v <- st_voronoi(x = st_union(mtq_c))
mtq_v <- st_intersection(st_cast(mtq_v), st_union(mtq))
mtq_v <- st_join(x = st_sf(mtq_v), y = mtq_c, join=st_intersects)
mtq_v <- st_cast(mtq_v, "MULTIPOLYGON")
plot(st_geometry(mtq_v), col='lightblue')
See file data_prep.R for data extraction.
# spatial data management
library(sf)
# cartography
library(cartography)
# smooth density computation
library(spatstat)
library(raster)
library(maptools)
# interactive mapping
library(mapview)
# load data
adm <- readRDS("data/iris_p13.rds")
feat_sir <- readRDS("data/sir_p13.rds")
feat_osm <- readRDS("data/osm_p13.rds")
nrow(feat_sir)
[1] 1066
nrow(feat_osm)
[1] 585
par(mar = c(0,0,1.2,0), mfrow = c(1,2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
plot(st_geometry(feat_sir), col = "red", pch = 20, add=T, cex = 0.5)
layoutLayer(title = "SIR", sources="", author="", scale = NULL, frame=FALSE)
legend(x = "topright", legend = "= 1 restaurant", pch = 20, pt.cex = 0.5, cex = 0.7, col = "red")
plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
plot(st_geometry(feat_osm), col = "red", pch = 20, add=T, cex = 0.5)
layoutLayer(title = "OSM", sources="", author="", scale = .5, north = T, frame=FALSE)
inter_osm <- st_intersects(adm, feat_osm)
inter_sir <- st_intersects(adm, feat_sir)
adm <- st_sf(adm[,1:2,drop = T],
n_osm = sapply(X = inter_osm, FUN = length),
n_sir = sapply(X = inter_sir, FUN = length),
geometry = st_geometry(adm))
par(mar = c(0,0,1.2,0), mfrow = c(1,2))
plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_sir", inches = 0.1)
layoutLayer(title = "SIR", sources="", author="", scale = NULL, frame=FALSE)
plot(st_geometry(adm), col = "ivory1", border = "ivory3",lwd =0.5,bg = "#FBEDDA")
propSymbolsLayer(adm, var = "n_osm", inches = 0.1)
layoutLayer(title = "OSM", sources="", author="", scale = .5, north = T, frame=FALSE)
# plot(adm$n_sir, adm$n_osm, asp = T, xlim = c(0,55), pch = 21, bg = "red")
# text(adm$n_sir, adm$n_osm, labels = adm$CODE_IRIS, cex = 0.5, pos = 2)
# x <-lm(formula = adm$n_osm~adm$n_sir)
# abline(a = x$coefficients )
Explore datasets interactively
The default is not bad:
mapview(feat_sir) + mapview(feat_osm, col.regions = "red")